#include "Astronomy/Precession.hpp"
#include "Astronomy/Transformation.hpp"


const PrecessionTable Precession::PrecessionTable1976
{
    {           0, 5038.7784, -1.07259, -0.001147, 0, 0},  // fi,ψ;psi;
    {   84381.448,         0, +0.05127, -0.007726, 0, 0},  // w,ω;omega;
    {           0,   +4.1976, +0.19447, -0.000179, 0, 0},  // P,
    {           0,  -46.8150, +0.05059, +0.000344, 0, 0},  // Q,
    {   84381.448,  -46.8150, -0.00059, +0.001813, 0, 0},  // E,ε;epsilon;
    {           0,  +10.5526, -2.38064, -0.001125, 0, 0},  // x;x;cai;
    {           0,   47.0028, -0.03301, +0.000057, 0, 0},  // pi(自导出)π;pi;
    {  629554.886, -869.8192, +0.03666, -0.001504, 0, 0},  // II(自导出)Π;Pi;
    {           0, 5029.0966, +1.11113, +0.000006, 0, 0},  // p;P;rho;
    {           0, 2004.3109, -0.42665, -0.041833, 0, 0},  // th,θ;theta;
    {           0, 2306.2181, +0.30188, +0.017998, 0, 0},  // Z,;Zeta;
    {           0, 2306.2181, +1.09468, +0.018203, 0, 0}   // z,ζ;zeta;

};

const PrecessionTable Precession::PrecessionTable2000
{

    {         0,  5038.478750,   -1.07259,  -0.001147,          0,          0},  // fi,ψ;psi;
    { 84381.448,    -0.025240,   +0.05127,  -0.007726,          0,          0},  // w,ω;omega;
    {         0,      +4.1976,   +0.19447,  -0.000179,          0,          0},  // P,
    {         0,     -46.8150,   +0.05059,  +0.000344,          0,          0},  // Q,
    { 84381.448,    -46.84024,   -0.00059,  +0.001813,          0,          0},  // E,ε;epsilon;
    {         0,     +10.5526,   -2.38064,  -0.001125,          0,          0},  // x;x;cai;
    {         0,      47.0028,   -0.03301,  +0.000057,          0,          0},  // pi(自导出)π;pi;
    {629554.886,    -869.8192,   +0.03666,  -0.001504,          0,          0},  // II(自导出)Π;Pi;
    {         0,   5028.79695,   +1.11113,  +0.000006,          0,          0},  // p;P;rho;
    {         0, 2004.1917476, -0.4269353, -0.0418251, -0.0000601, -0.0000001},  // th,θ;theta;
    {+2.5976176, 2306.0809506, +0.3019015, +0.0179663, -0.0000327, -0.0000002},  // Z,;Zeta;
    {-2.5976176, 2306.0803226, +1.0947790, +0.0182273, +0.0000470, -0.0000003}   // z,ζ;zeta;

};

const  PrecessionTable Precession::PrecessionTableP03
{
    {           0, 5038.481507,  -1.0790069,  -0.00114045, +0.000132851, -9.51e-8 },   // fi,ψ;psi;
    {84381.406000,   -0.025754,  +0.0512623,  -0.00772503,     -4.67e-7, +3.337e-7},   // w,ω;omega;
    {           0,    4.199094,  +0.1939873,  -0.00022466,     -9.12e-7, +1.20e-8 },   // P,
    {           0,  -46.811015,  +0.0510283,  +0.00052413,     -6.46e-7, -1.72e-8 },   // Q,
    {84381.406000,  -46.836769,  -0.0001831,  +0.00200340,     -5.76e-7, -4.34e-8 },   // E,ε;epsilon;
    {           0,   10.556403,  -2.3814292,  -0.00121197, +0.000170663, -5.60e-8 },   // x;x;cai;
    {           0,   46.998973,  -0.0334926,  -0.00012559,     +1.13e-7, -2.2e-9  },   // pi(自导出)π;pi;
    { 629546.7936, -867.957580,  +0.1579920,  -0.00053710,  -0.00004797, +7.2e-8  },   // II(自导出)Π;Pi;
    {           0, 5028.796195,  +1.1054348,  +0.00007964, -0.000023857, +3.83e-8 },   // p;P;rho;
    {           0, 2004.191903,  -0.4294934,  -0.04182264,    -7.089e-6, -1.274e-7},   // th,θ;theta;
    {    2.650545, 2306.083227,  +0.2988499,  +0.01801828,    -5.971e-6, -3.173e-7},   // Z,;Zeta;
    {   -2.650545, 2306.077181,  +1.0927348,  +0.01826837, -0.000028596, -2.904e-7}    // z,ζ;zeta;
};

double Precession::CalcEntry(double T, const double entry[], const size_t entrySize) noexcept
{
    double result{ 0 };

    for (size_t i = 0; i < entrySize; i++, T *= T)
    {
        result += entry[i] * T;
    }
    return Transformation::ArcsecToRadians(result);
}

std::tuple<const PrecessionTable*, size_t> Precession::GetPrecessionTablePtr(const PrecessionType type) noexcept
{
    switch (type)
    {
    case PrecessionType::IAU1978:
        return  { &PrecessionTable1976 ,4 };
    case PrecessionType::IAU2000:
        return { &PrecessionTable2000,6 };
    //case PrecessionType::IAUP03:
    //    return { &PrecessionTableP03,6 };
    default:
        break;
    }
    return { &PrecessionTableP03,6 };
}

Coordinate3D Precession::CalcEquatorConversion(const Coordinate3D & llr, double phi, double sigma, double zeta) noexcept
{
    const double Alpha{ llr.X };
    const double Delta{ llr.Y };

    //var cosW = cos(llr[1]), cosH = cos(th);
    //var sinW = sin(llr[1]), sinH = sin(th);
    const double cosDelta = cos(Delta);
    const double sinDelta = sin(Delta);

    const double cosphi = cos(phi);
    const double sinphi = sin(phi);

    // var Z = prece(t,'Z', mx) + llr[0];
    // cos(Z)
    const double cosAlphaplussigma = cos(Alpha + sigma);

    //var A = cosW * sin(Z);
    const double A = cosDelta * sin(sigma);

    //var B = cosH * cosW * cos(Z) - sinH * sinW;
    const double B = cosphi * cosDelta * cosAlphaplussigma - sinphi * sinDelta;

    //var C = sinH * cosW * cos(Z) + cosH * sinW;
    const double C = sinphi * cosDelta * cosAlphaplussigma + cosphi * sinDelta;

    //return new Array(rad2mrad(atan2(A, B) + z), asin(C), llr[2]);
    return{ Transformation::MapTo0To2PIRange(atan2(A, B) + zeta),asin(C),llr.Z };
}

Coordinate3D Precession::CalcEclipticConversion(const Coordinate3D &llr, double eta, double pi, double p) noexcept
{
    const double Lambda{ llr.X };
    const double Beta{ llr.Y };

    const double cosBeta = cos(Beta);
    const double sinBeta = sin(Beta);

    const double coseta = cos(eta);
    const double sineta = sin(eta);

    const double sinpiminusLambda = sin(pi - Lambda);

    const double A = (coseta * cosBeta * sinpiminusLambda) - (sineta * sinBeta);
    const double B = cosBeta * cos(pi - Lambda);
    const double C = (coseta * sinBeta) + (sineta * cosBeta * sinpiminusLambda);

    return{ Transformation::MapTo0To2PIRange(p + pi - atan2(A, B)),asin(C),llr.Z };
}


Coordinate3D Precession::ToDateEquator(const double T, const Coordinate3D llr, const PrecessionType type) noexcept
{
    auto [pPrecessionTable, Count] = GetPrecessionTablePtr(type);

    double phi   = CalcEntry(T, pPrecessionTable->Phi, Count);
    double sigma = CalcEntry(T, pPrecessionTable->Sigma, Count);
    double zeta  = CalcEntry(T, pPrecessionTable->Zeta, Count);

    return  CalcEquatorConversion(llr, phi, sigma, zeta);
}

Coordinate3D Precession::ToJ2000Equator(const double T, const Coordinate3D llr, const PrecessionType type) noexcept
{
    auto [pPrecessionTable, Count] = GetPrecessionTablePtr(type);

    double sigma = -CalcEntry(T, pPrecessionTable->Sigma, Count);
    double phi   = -CalcEntry(T, pPrecessionTable->Phi, Count);
    double zeta  = -CalcEntry(T, pPrecessionTable->Zeta, Count);

    return  CalcEquatorConversion(llr, sigma, phi, zeta);
}

Coordinate3D  Precession::ToDateEcliptic(const double T, const Coordinate3D llr, const PrecessionType type) noexcept
{
    auto [pPrecessionTable, Count] = GetPrecessionTablePtr(type);


    double epsilon = CalcEntry(T, pPrecessionTable->Epsilon, Count);//E
    double lambda  = CalcEntry(T, pPrecessionTable->Lambda, Count);//x;
    double cmega   = CalcEntry(T, pPrecessionTable->Omega, Count);//w;
    double psi     = CalcEntry(T, pPrecessionTable->Psi, Count);//fi;

    auto result{ llr };

    result.X += psi;
    result.EclipticToEquatorial(cmega);
    result.X += -lambda;
    result.EclipticToEquatorial(-epsilon);

    return result;
}

Coordinate3D  Precession::ToJ2000Ecliptic(const double T, const Coordinate3D llr, const PrecessionType type) noexcept
{
    auto [pPrecessionTable, Count] = GetPrecessionTablePtr(type);


    double epsilon = CalcEntry(T, pPrecessionTable->Epsilon, Count);//E
    double lambda  = CalcEntry(T, pPrecessionTable->Lambda, Count);//x;
    double cmega   = CalcEntry(T, pPrecessionTable->Omega, Count);//w;
    double psi     = CalcEntry(T, pPrecessionTable->Psi, Count);//fi;

    auto result{ llr };

    result.EclipticToEquatorial(epsilon);
    result.X += lambda;
    result.EclipticToEquatorial(-cmega);
    result.X -= psi;

    return result;
}



